Structural Origin of Dynamical Heterogeneity in a Binary Lennard- Jones Liquid 
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We use isoconfigurational ensemble simulations to analyze dynamical heterogeneity (DH) and 
structural heterogeneity (SH) as a function of temperature T in the supercooled Kob-Andersen 
binary Lennard- Jones liquid. At high T we identify SH that is larger than the observed DH. As T 
decreases, the DH grows but remains smaller than or comparable in size to the SH at all T studied. 
Our results thus elucidate the structural origin of DH in this important model glass-former. 
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Glass-forming liquids are remarkable for the extraordi- 
nary sensitivity of their transport properties to changes 
in state variables, such as temperature T [T]. Proper- 
ties such as viscosity are commonly found to vary over 
14 orders of magnitude in supercooled liquids between 
the melting temperature and the glass transition. Yet in 
the same interval of T, it is also typical that only mod- 
est changes occur in the average liquid structure. One 
of the central questions in the study of the glass transi- 
tion is whether this enormous dynamical response can be 
understood in terms of structural change [5] . 

Much recent interest has focussed on the emergence 
and growth of dynamical heterogeneity (DH) in super- 
cooled liquids, that is, spatially extended domains in 
which molecules are more or less mobile, relative to the 
bulk average |3]. The growth of these dynamical do- 
mains as T decreases seems to occur in the absence of 
a growing structural length scale. At the same time, 
the simulation work of Widmer-Cooper, Harrowell and 
coworkers has shown that key aspects of DH are indeed 
structural in origin [11|5J[3]. They do so through the use 
of the "isoconfigurational (IC) ensemble", a simulation 
procedure in which a given liquid configuration is ana- 
lyzed by conducting a set of runs all initiated from the 
same configuration, but in which particle velocities are 
randomized [4^. Analysis of the "dynamic propensity", a 
particle's displacement averaged over all the runs of the 
IC ensemble, reveals spatial heterogeneity that can only 
be due to structural properties of the initial configura- 
tion, since the effect of the velocity initial conditions has 
been averaged out. 

In this Letter, we evaluate the dynamic propensity 
in the Kob-Andersen binary Lennard- Jones liquid [7], a 
model system that has received much attention in simu- 
lations of glass- forming liquids in general [S] , and in the 
analysis of DH in particular [9l [TOl [TTJ [12]. We also an- 
alyze the structure of the liquid via IC averaging, and 
show that structural heterogeneity (SH) exists at high 
T that is greater in size than the DH occurring at the 
same T. Further, as T decreases and the DH grows, its 
size remains smaller than, or comparable to, the SH that 
we find at all T studied here. Although there has been 
significant recent progress in identifying local structural 
properties (e.g. potential energy, soft vibrational modes) 
correlated to DH [HI HSl HH HI] , to our knowledge there 



are no reports of SH that exceeds the size of DH in the 
glass- forming regime, as we show here. 

Our model system is the Kob-Andersen liquid, con- 
sisting of an 80:20 mixture of = 8788 A and B parti- 
cles, interacting via a potential Vap = 4,ea0[{crap/i")^'^ — 
(fQ/s/?")^] with g {A,B}. All quantities are reported 
here in reduced units, with length, energy and time given 
relative to aAA, ^aa and (wcr^^/48eyi^)^/^ respectively, 
where m is the mass of the particles. The potential 
parameters are aAA = 1-0, ^aa = 1-0, ctbb = 0.88, 
Ebb = 0.5, aAB = 0.8 and eAB = 1-5 [7J. The poten- 
tial is truncated and shifted at a cutoff radius of 2.5. All 
simulations are conducted in a cubic cell with periodic 
boundary conditions at a fixed volume V to give a den- 
sity of 1.2. The simulation time step is At = 0.01. In all 
cases below, we restrict our attention to the properties 
of the A particles. 

We study the liquid using the IC ensemble method at 
T = 1.0, 0.6, 0.5 and 0.466. At each T we generate 10 
independent starting configurations. We first equilibrate 
a random configuration at T = 5.0 for at least 28 000 time 
steps, then reset T to the desired value, controlling T 
throughout with a Berendsen thermostat. Each system is 
equilibrated for at least 20Ta where Tq, is the value of the 
a-relaxation time at that T. We then use each starting 
configuration to initiate M = 500 runs of an IC ensemble 
by randomizing the velocities of all particles according 
to a Maxwell-Boltzmann distribution, while leaving the 
particle coordinates unchanged. The IC ensemble runs 
are carried out in the microcanonical ensemble. 

Let r'^{i,k,t) be the squared displacement of the i-th 
A particle at time t in run k of an IC ensemble. The dy- 
namic propensity of each A particle at time t is the value 
of (r2)i, = M-iX;fli^^(^fc,0- We evaluate (rf);^ for 
each A particle as a function of t, and quantify the spa- 
tial correlations of this quantity, in terms of mobile and 
immobile domains, via a cluster analysis [TH [TS]. We 
identify the particles having the highest 10% of (rf)ic 
values at a given t, and then find the clusters formed by 
this subset. Clusters are defined by the criterion that 
two particles of the subset that are also within r = 1.4 
of one another (the position of the first minimum of the 
A-A radial distribution function) in the initial configu- 
ration are assigned to the same cluster. The number- 
averaged mean cluster size of a set of Nc clusters is 
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FIG. 1: Mean cluster size as a function of t for (a) mobile 
clusters, (b) immobile clusters, (c) clusters of A particles with 
low B coordination, and (d) clusters of A particles with high 
B coordination, all evaluated in the IC ensemble. 



S — (l/Nc) J^n nC{n), where C(n) is the number of clus- 
ters of size n. We evaluate S for the clusters of mobile 
particles defined above, and denote it as Sm- We con- 
duct the same analysis on the lowest 10% of (rf )ic values, 
and find the mean cluster size of this immobile subset Sj. 
Figs, [ija) and (b) show the t dependence of Sm and Sj, 
where the data are averaged over the 10 starting config- 
urations used at each T. 

For comparison, we also evaluate S for the DH that 
occurs in single simulation runs. That is, we separately 
analyze each of the M = 500 runs of one IC ensemble at 
each T, and evaluate, as a function of t, the mean cluster 
size Sm of the clusters formed by the particles having the 
largest 10% of displacements as measured from their po- 
sition in the starting configuration. These mobile clusters 
correspond to the "strings" documented e.g. in Ref. [H]. 
We then obtain the average of the Sm curves over all 500 
runs [Fig. [2]^ a)]. The corresponding cluster-size analysis 
of the smallest 10% of displacements in individual runs 
gives Si as a function of t [Fig. [2j^b)]. 




FIG. 2: Mean cluster size as a function of t for (a) mobile 
clusters and (b) immobile clusters, as measured in single sim- 
ulations runs. 



Fig. [T];a,b) and Fig. [2] allow us to compare the DH 
as revealed by both IC and conventional averaging, for 
both mobile and immobile domains. The t dependence 
of all curves follows the behavior for DH found in earlier 
work (see e.g. Ref. [IB|). At small t, S has the value 
expected for a random choice of 10% of the A particles 
(approximately So — 2.135), consistent with no spatial 
correlations. However, on the time scale of structural re- 
laxation a maximum occurs, indicating significant clus- 
tering of mobile and immobile particles. At large t, the 
DH begins to dissipate and S decreases toward So- 

The monotonic increase in S"^^^ (the maximum value 
of S) as T decreases quantifies the growth of DH on cool- 
ing. We denote the maximum value of Sm in Fig.[lja) as 
S]^^'^; similarly, S'™^'^, S'^^^'^ and 5,'°^'' denote the maxima 
in Figs. [l](b), [2](a) and |2jb) respectively. The T depen- 
dence of 5'™^'^ for all data in Figs. [l]and|2]is plotted in 
Fig. |3] as a function of (T - Te)/Tc, where = 0.435 
is the critical temperature of mode coupling theory for 
the Kob- Andersen liquid. Fig. [3] shows that the sizes of 
the mobile and immobile clusters found using the IC en- 
semble are as much as an order of magnitude larger than 
the DH found when analyzing single runs. Also, the T 
dependence of S^^^ is quite different for the two kinds 
of averaging. 5,™^'^ and Sf^'^^ follow a power law f9| IT?], 
while S'^l^^ and SY^^ do not. Indeed, the most notable 
behavior in Fig. [3] is that Sft"^ and Sf^^ both initially 
grow faster than a power law on cooling, but then at the 
lowest T their growth slows considerably. 

To explore this behavior, in Fig. |4] (top panels) we vi- 
sualize the mobile and immobile domains of the dynamic 
propensity for one starting configuration at each T, at the 
value of t corresponding to S^f^, where the clusters are 
most prominent [14]. Particles in the top (bottom) 50% 
of {r1)ic values are represented as green (red) spheres, 
with each sphere plotted at the position of the particle 
in the initial configuration. The radius of each sphere 
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(T-T,)/T, 

FIG. 3: Maximum value 5'™'"' of the mean cluster size for 
each of the curves plotted in Figs, [l] and [2] as a function of 
(T—Tc)/T, where Tc = 0.435. Shown are the maximum values 
of Si (filled squares), Sm (filled circles), Si (open squares), 
Sm (open cicles), Sh (up triangles) and Sl (down triangles). 
So ~ 2.135 has been subtracted from all the data to refiect 
the value of 5'™'''' in excess of its value for a random choice of 
10% of the A particles. For the open symbols, the standard 
deviation of the data over the 10 starting configurations is 
less than or comparable to the symbol size. 

represents the rank order of {rf)ic- the larger a green 
(red) sphere is, the larger (smaller) is its value of {rf)ic- 
(See caption of Fig. |4]for complete details.) These visu- 
alizations illustrate the increasing spatial correlation in 
the dynamic propensity as T decreases. At T = 1.0 the 
arrangement of red and green spheres is nearly random, 
while at the lowest T, very large mobile and immobile 
domains have emerged. At both T — 0.5 and 0.466, the 
domains are strikingly large, and are comparable to the 
system size. This suggests that finite- size effects may be 
responsible for the low T behavior of SfY'^ and S™^^ in 
Fig. |3] We return to this issue below. 

We next test for the presence of SH in the liquid, by 
analyzing a structural property of the system using the 
same IC averaging employed to find the dynamic propen- 
sity |14j . Specifically, we examine the chemical compo- 
sition of the nearest neighbor environment of the A par- 
ticles. Let n(i,fc,t) be the number of B particles found 
within a distance of r = 1.2 (the first minimum of the 
A-B radial distribution function) of the z-th A particle 
at time t in run k of an IC ensemble. We define the "B- 
coordination propensity" of each A particle at time t as 
{ni)ic = M~^^^^^^n{i,k,t). We show in Fig. Il|c) the 
time dependence of the mean cluster size Sl of the 10% 
of A particles with the lowest values of {ni)ic- Fig. ^d) 
shows Sh, the mean cluster size of the 10% of A particles 
with the highest values of (n,i)ic. 

The time evolution of Sl and 5*^^ in Fig. [T|is similar to 
Sm and Sj with one surprising exception: The values of 
^max g^^^ S'^^^ (shown in Fig. |3]) are nearly independent 



of T, and are larger than or comparable to 5'^'^'^ and Sf^^^ 
for all T. The spatial variation of (ni)ic is visualized in 
the bottom panels of Fig. |4]for the same initial configura- 
tions shown in the top panels, with t chosen at the time 
of S'^'^'^. These visualizations confirm that the size of the 
domains with low and high B coordination remain large 
at all T, including at high T where the mobile and im- 
mobile domains are an order of magnitude smaller |18j . 
Further, the locations of the mobile and immobile do- 
mains that emerge on cooling approximately correspond 
with the domains of low and high B coordination that 
are prominent at all T. This spatial correspondence is 
consistent with expectation: an A particle with low B 
coordination will be less tightly bound by it's neighbors, 
and thus potentially more mobile, than one with high B 
coordination. 

Hence, averaging a local structural property (here, the 
B coordination of A particles) in the IC ensemble reveals 
SH the size of which is insensitive to T over the range 
studied here. Of itself, this behavior is entirely consis- 
tent with the established view that structural correlation 
lengths do not grow significantly in most supercooled liq- 
uids. However, what is new in our results is that the ob- 
served SH is larger than or comparable to the DH by any 
of our measures (i.e. via IC or conventional averaging). 
In general, the increasing sensitivity of local dynamics to 
local structure as T decreases makes it inevitable that 
DH will occur in any liquid on a scale at least up to the 
size of any SH that is present. This scenario, in which DH 
is a dynamical consequence of an underlying SH, appears 
to apply in the range of T studied here. 

We emphasize that since our main results are based 
on isoconfigurational averaging, the structural triggers 
for individual correlated dynamical events occurring in 
a single simulation run (e.g. the "strings" of Ref. [3]) 
are not specifically addressed here. Our results identify 
the structural properties that make such events more or 
less probable, and which induce the mobile and immobile 
domains visualized in Fig. [4] 

An open question raised by our results is the role of 
finite size effects. The size of the domains observed in 
Fig. |4] certainly suggests that larger systems are required 
to fully understand the nature of the observed correla- 
tions, and their potential impact on other observables, 
both static and dynamic. For example, it is possible that 
in larger simulations the size of the DH may grow be- 
yond that of the SH for T < 0.5. This would indicate 
an interesting crossover of the mechanism for DH, from 
one controlled by SH, to one allowing growth beyond the 
scale set by SH. However, our conclusions as they apply 
to DH for T > 0.5 are unlikely to change in simulations 
of larger systems, since the DH becomes smaller than the 
system size in this regime. 
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FIG. 4: Spatial variation of (rf )ic (top panels) and {ni)ic (bottom panels) at each T. To make each top panel, the values of 
(r^)ic, evaluated at the time of the maximum of Sm, are assigned to each A particle at its position in the initial configuration 
of the IC ensemble. These values are sorted and assigned an integer rank Ri from 1 to A'^, from smallest to largest. Each A 
particle is then plotted as a green sphere of radius a = Rmin exp{ [{Ri — N)/{1 — N)] log(_Rmax/ Rmin)}, where i?max = 0.5 and 
Rmin = 0.01. The ranks Ri are then reversed (i.e. assigned from largest to smallest), and each A particle is also plotted as 
a red sphere of radius a. The color observed for each particle therefore indicates which of the green or red spheres is larger. 
The result presents the rank of {r^)ic on an exponential scale, such that the largest green spheres represent the most mobile A 
particles, and the largest red spheres the most immobile. The bottom panels are created in exactly the same way as the top 
panels, but with {r^)ic replaced by {ni)ic, and where the time is chosen to be the maximum of Sl at each T. In the bottom 
panels, the largest yellow spheres represent the A particles with the lowest B coordination, and the largest blue spheres the A 
particles with the highest B coordination. 
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